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We analyze the realization of a quantum-walk search algorithm in a passive, linear optical network. 
The specific model enables us to consider the effect of realistic sources of noise and losses on the 
search efficiency. Photon loss uniform in all directions is shown to lead to the rescaling of search time. 
Deviation from directional uniformity leads to the enhancement of the search efficiency compared 
to uniform loss with the same average. In certain cases even increasing loss in some of the directions 
can improve search efficiency. We show that while we approach the classical limit of the general 
search algorithm by introducing random phase fluctuations, its utility for searching is lost. Using 
numerical methods, we found that for static phase errors the averaged search efficiency displays a 
damped oscillatory behaviour that asymptotically tends to a non-zero value. 



I. INTRODUCTION 



The generalization of random walks for quantum sys- 
tems QJ proved to be a fruitful concept [2| attracting 
much recent interest. Algorithmic application for quan- 
tum information processing is an especially promising 
area of utilization of quantum random walks (QRW) [||. 

In his pioneering paper Q Grover presented a quan- 
tum algorithm that can be used to search an unsorted 
database quadratically faster than the existing classical 
algorithms. Shenvi, Kempe and Whaley (SKW) [j| pro- 
posed a search algorithm based on quantum random walk 
on a hypcrcubc, which has similar scaling properties as 
the Grover search. In the SKW algorithm the oracle is 
used to modify the quantum coin at the marked vertex. 
In contrast to the Grover search, this algorithm generally 
has to be repeated several times to produce a result, but 
this merely adds a fixed overhead independent of the size 
of the search space. 

There are various suggestions and some experiments 
how to realize quantum walks in a laboratory. The 
schemes proposed specifically for the implementation of 
QRWs include ion traps @, nuclear magnetic resonance 
pj (also experimentally verified @), cavity quantum elec- 
trodynamics [{J fTol | , o ptical lattices 11 1, optical traps 
[HI, optical cavity [IH, and classical optics More- 
over, the application of standard general logic networks 
to the task is always at hand [H, [l(| . 

The idea of the scattering quantum random walk 
(SQRW) [13 was proposed as an answer to the question 



that can be posed as: how to realize a coined walk by a 
quantum optical network built from passive, linear opti- 
cal elements such as beam splitters and phase shifters? It 
turned out that such a realization is possible and, in fact, 
it leads to a natural generalization of the coined walk, the 
scattering quantum random walk [la ]. The SQRW on the 
hypercube allows for a quantum optical implementation 
of the SKW search algorithm Q . Having a proposal for 
a physical realization at hand we are in the position to 
analyze in some detail the effects hindering its successful 
operation. 



Noise and decoherence strongly influence quantum 
walks. For a recent review on this topic see P~9[- The 
first investigations in this direction indicated that a small 
amount of decoherence can actually enhance the mixing 
property [2(j. For a continuous QRW on a hypercube 
there is a threshold for decoherence, beyond which the 
walk behaves classically [2l| . Kosik et al analyzed SQRW 
with randomized phase noise on a d dimensional lattice 
[22| . The quantum walk on the line has been studied 
by several authors in the linear optical context, with the 
emphasis on the effect of various initial states, as well 
as on the impact of decoherence [2^, [24[ . The quantum 
random walk search with imperfect gates was discussed 
in some detail by Li et al [25j, who have considered the 
case when the Grover operator applied in the search is 
systematically modified. Such an imperfection decreases 
the search probability and also shifts its first maximum 
in time. 

In this paper we analyze the impact of noise on the 
SKW algorithm typical for the experimental situations 
of the SQRW. In particular, first we focus on photon 
losses and show that, somewhat contradicting the naive 
expectation, non-trivial effects such as the enhancement 
of the search efficiency can be observed. As a second type 
of errors we study randomly distributed phase errors in 
two complementary regimes. The first regime is charac- 
terized by rapid fluctuation of the optical path lengths, 
that leads to the randomization of phases for each run 
of the algorithm. We show that the classical limit of the 
SKW algorithm, reached by increasing the variance of 
the phase fluctuations, does not correspond to a search 
algorithm. In the other regime, the stability of the opti- 
cal path lengths is maintained over the duration of one 
run, thus the errors are caused by static random phases. 
This latter case has not yet been considered in the con- 
text of QRWs. We found that static phase errors bring 
a significantly different behaviour compared to the case 
of phase fluctuations. Under static phase errors the al- 
gorithm retains its utility, with the average success prob- 
ability displaying a damped oscillatory behaviour that 
asymptotically tends to a non-zero constant value. 

The paper is organized as follows. In the next section 
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we introduce the scattering quantum walk search algo- 
rithm. In section IIIH we derive analytic results for the 
success probability of search for the case when a single 
coefficient describes photon losses independent of the di- 
rection. In section IIVI we turn to direction dependent 
losses, and present estimations of the success probability 
based on analytical calculations and numerical evidence. 
In section [V] phase noise is considered and consequences 
for the success probability are worked out. Finally, we 
conclude in Sec. [VI] 



II. THE SCATTERING QUANTUM WALK 
SEARCH ALGORITHM 

The quantum walk search algorithm is based on the 
generalized notion of coined quantum random walk 
(CQRW), allowing the coin operator to be non- uniform 
accross the vertices. In the early literature the coin is 
considered as position (vertex) independent. The CQRW 
is defined on the product Hilbert space H. = 7i c <S> H G , 
where 7i c refers to the quantum coin, and TL repre- 
sents the graph on which the walker moves. The discrete 
time-evolution of the system is governed by the unitary 
operator 

U = SC, (1) 

where C is the coin operator which corresponds to flip- 
ping the quantum coin, and S is the step or translation 
operator that moves the walker one step along some out- 
going edge, depending on the coin state. Adopting a 
binary string representation of the vertices V of the un- 
derlying graph G = (V,E), the step operator S (a per- 
mutation operator of the entire Hilbert space TL) can be 
expressed as 

ri-1 

S=J2J2\d,x®e dx )(d,x\ . (2) 
d=o xev 

In ^ x denotes the vertex index. Here, and in the rest 
of this paper we identify the vertices with their indices 
and understand V as the set of vertex indices. The most 
remarkable fact about S is that it contains all information 
about the topology of the graph. In particular, the actual 
binary string values of e dx are determined by the set of 
edges E. This is accomplished by the introduction of 
direction indices d, which run from to n—1 in case of the 
n regular graphs which are used in the search algorithm. 

To implement the scattering quantum random walk 
on an n regular graph of N nodes, identical ri-multiports 
[HI, H3| are arranged in columns each containing N mul- 
tiports. The columns are enumerated from left to right, 
and each row is assigned a number sequentially. The ini- 
tial state enters on the input ports of multiports in the 
leftmost column. The output and input ports of multi- 
ports of neighbouring columns j and j+ are then indexed 
suitably and connected according to the graph G. 



For the formal description of quantum walks on arrays 
of multiports, we propose to label every mode by the row 
index and input port index of its destination multiport. 
We note that an equally good labelling can be defined 
using the row index and output port index of the source 
multiport. To describe single excitation states, we use the 
notation \d, x) where the input port index of the destina- 
tion multiport is d = 0, 1, . . . , n ,— I, and the row index is 
x = 0, 1, . . . , N — 1. Thus the total Hilbert space can ef- 
fectively be separated into some product space TL C ®TL G . 
To be precise, the additional label j would be necessary 
to identify in which column the multiport is, however, 
we think of the column index as a discrete time index, 
and drop it as an explicit label of modes. Thus a time- 
evolution U = SC can be generated by the propagation 
through columns of multiports. 

A quantum walk can be realized in terms of the basis 
defined using the destination indices, and we shall term 
it "standard basis" through this section. First, we re- 
call that an n-multiport can be fully characterized by an 
SU(n) transformation matrix C. The effect of such mul- 
tiport on single excitation states G TL° is given by 
the formula, 

n— 1 n — 1 

|V) =5> d |d) -+ C dka k \d), (3) 

d=0 d,k=0 

where \d) denotes the single photon state with the pho- 
ton being in the d mode, i.e. \d) = |0) . . . |l) d . . . |0) n _ 1 . 
We note, that a multiport with any particular transfor- 
mation matrix C can be realized in a laboratory [28| . To 
simplify calculations it may be beneficial to choose an 
indexing of input and output ports such that the con- 
nections required to realize the graph G can be made in 
such way that each input port has the same index as the 
corresponding source output port. Therefore the label 
d can stay unique during "propagation." We emphasize 
that this is not a necessary assumption for a proper def- 
inition of SQRW, but an important property that makes 
also easier to see that SQRWs are a superset of general- 
ized CQRWs. This indexing of input and output ports 
for walks on a hypercube is depicted on Fig. [T^,, with 
some of the actual connections illustrated for a (three 
dimensional) cube on Fig. [TJd. 

Considering an array of identical multiports, an ar- 
bitrary input state undergoes the transformation by the 
same matrix C for every x. Let the output port d of mul- 
tiport x be connected to multiport x © e dx in the next 
row. Thus the mode labelled by the source indices d and 
x, is labelled by d and x®e dx in terms of the destination 
indices. Therefore, effect of propagation in terms of our 
standard basis is written, 

^ a dx \d, x) -> ^ C dk a k x |o, X © 6dx )■ (4) 

d,x dkx 

Comparing this formula with Eqs. © and we see that 
this formula corresponds to a U = SC = 5(Co<E>l) trans- 
formation where Co is generated by the matrix C. Due 
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FIG. 1: a) Illustration of the labelling of input and output 
ports of multiports used for the realization of the walk on 
the n dimensional hypercube. b) Schematic depiction of the 
setup of the SQRW implementation of the SKW algorithm 
for n = 3, with the marked node being x t = 001. For clarity, 
a unique line pattern is associated with each input (output) 
port index. 



to the local nature of the realization of the coin opera- 
tion, it is straight-forward to realize position dependent 
coin operations, such as the one required for the quantum 
walk search algorithm. 

In particular, the SKW algorithm [f| is based on the 
application of two distinct coin operators, e.g. 



Co — G, 

Ci = -i, 



(5a) 
(5b) 



where G is the Grover inversion or diffusion o per ator 
G := -t + 2\s c ){s c \, with \s°) = 1/V^ELi \ d ) & In 
the algorithm, the application of the two coin operators 
is conditioned on the result of oracle operator O. The 
oracle marks one Xt as target, hence the coin operator 
becomes conditioned on the node: 



C = C ®1 + (C 1 -C )®\x t )(x t \ 



(6) 



When n is large, the operator U' := SC can be regarded 
as a perturbed variation of U = S(Co ® 1). The con- 
ditional transformation §Q is straight-forward to imple- 
ment in the multiport network. For the two coins © 
one has to use a simple phase shifter at position x t in 
every column of the array, and a multiport realizing the 
Grover matrix G at every other position. The connection 
topology required to implement a walk on the hyper- 
cube is such that in the binary representation we have 
e c i = ... 1 ... with 1 being at the d'th position, i.e. 
e c i = 2 d . See Fig. [T|d for a schematic example, when 
x t = 001. 

The above described scheme to realize quantum walks 
in an array of multiports using as many columns as the 
number of iterations of U can be reduced to only a single 
column. To do this, one simply needs to connect the 
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FIG. 2: Schematic illustration of the photon loss model being 
used. The losses suffered by each output mode are repre- 
sented by fictitious beam-splitters with transmittivities r/d- 
The beam-splitters incorporate the combined effect of imper- 
fections of the multiport devices, and effects influencing the 
state during propagation between the multiports (e.g. scat- 
tering and absorption). 



output ports back to the appropriate input ports of the 
destination multiport in the same column. This feed- 
back setup is similar to the one introduced in Ref . [IU . 



III. UNIFORM DECAY 

We begin our analysis of the effect of errors on the 
quantum walk search algorithm by concentrating on pho- 
ton losses. In an optical network, photon losses are usu- 
ally present due to imperfect optical elements. An effi- 
cient model for linear loss is to introduce fictitious beam- 
splitters with transmittances corresponding to the effec- 
tive transmission rate (see Fig. [5]). 

The simplest case is when all arms of the multiports 
are characterized by the same linear loss rate rj. The op- 
erator describing the effect of decay on a single excitation 
density operator can then be expressed as 



ri 2 g+{l-i] 2 ) 



(7) 



The total evolution of the system after one iteration may 
be written as g — > T>{U gU^). It is important to note that 
with the introduction of this error, the original Hilbcrt 
space H G of one-photon excitations must be extended 
by the addition of the vacuum state |0). The action of 
the SQRW evolution operator U on the extended Hilbcrt 
space follows from the property U\0) = |0). Due to the 
nature of Eq. ([7]) and the extension of U, one can see that 
the order of applying the unitary time step and the error 
operator D can be interchanged. Therefore, over t steps 
the state of the system undergoes the transformation 

q -> rftfoU* + (1 - v 2t ) |0> <0| = V^U'gU^). (8) 

To simplify calculations, we introduce a linear (but 
non-unitary) operator to denote the effect of the noise 
operator T> on the search Hilbert space: 



D\i>)=r,\ip). 



(9) 



This operator is simply a multiplication with a number. 
It is obviously linear, however, for 77 < 1 not unitary. 
The operator D does not describe any coherence damping 
within the onc-photon subspace, since it only uniformly 
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decreases the amplitude of the computational states and 
introduces the vacuum. Since all final statistics are gath- 
ered from the search Hilbcrt space Tt G , it is possible to 
drop the vacuum from all calculations, and incorporate 
all information related to it into the norm of the remain- 
ing state. In other words, we can think of DU as the 
time step operator, and relax the requirement of normal- 
ization. Using this notation, the effect of t steps is very 
straight-forward to express: 



rfV l \ 



(10) 



This formula indicates that inclusion of the effect of uni- 
form loss may be postponed until just before the final 
measurement. The losses, therefore, may simply be in- 
cluded in the detector efficiency (using an exponential 
function of the number of iterations). 

Applying the above model of decay to the quantum 
walk search algorithm we define the new step operator 
U" = DU' , and write the final state of the system after 
t steps as 

(u"Y l^o) = 

if cos(^t) |Vo) - rf sin(^i) + rfO \-j=J \f) . 

(11) 

Adopting the notation of Ref. Q , the probability of mea- 
suring the target state |x = 0) at the output after t steps 
can be expressed as 



PnM) = Od,0|(C/")Vo>r 

= v 2t sin 2 (u' t) | (R, 0| i/fi) | 2 + 2"V' cos 2 (uj' Q t) 

+ 0(1/2"). (12) 

We know from Ref. [5] that | (i?, 0| t/>i) | 2 = l/2-0(l/n). 
Since an overall exponential drop of the success proba- 
bility is expected due to the r] 2t factor, we search for the 
maximum tf be before the ideal time-point \u/' \t = tt/2. 
This guarantees that sin 2 (wgt/) is finite, therefore due to 
the 2~ n factor for large n the second term can be omitted, 
and it is sufficient to maximize the function 



p n ( V , t) = ?? 2t sin 2 K0 (1/2 - 0(l/n)) 



(13) 



with respect to t. After substituting the result \u>' \ = 
l/^/2 n - l [l - 0(1/ n)] ± 0(n 3 / 2 /2") from Ref. 0], these 
considerations yield the global maximum at tf = 

During operation 



V2™ 31 
we set 



acot(- \nr]\ / 2 n - 1 ) + 0(l/n) 



To simplify the upcoming formulae, we introduce the 
variables 



- In 77 \P2F- 
log 2 (l - 77). 



(15) 
(16) 



The variable e can be regarded as a logarithmic trans- 
mission parameter (the ideal case corresponds to e = 00, 
and complete loss to e = 0). When e is sufficiently large, 
the expression — In 77 can be approximated to first order 
in 2~ e and we obtain 



-e+n/2^1/2 



(17) 



Upon substituting t m into (|13p we can use the new vari- 
able x to express the sine term as 

sin 2 (wot m ) = sin 2 (|wQ|7: m ) = sin 2 [acota;(l + 0(l/n))] = 

1 2x acot x . . . acot 2 x . , 9. 

0{l/n) + ——^0{l/n 2 ). (18) 



1 + x 



1 + x 2 



1 + x 2 



Thus for the maximum success probability p™ ax (?7) = 
Pn(Vitm) we obtain 



max 



in) 



1 



\-0{l/n) 



x acot x 0(l/n) 



(19) 

This formula is our main result for the case of uniform 
photon losses. In the large n limit it gives the approx- 
imate performance of the SKW search algorithm as a 
function of the transmission rate and the size of the 
search space. Since x acot x is bounded in x, the accu- 
racy of the term in brackets is bounded by 0(l/n). The 
most notable consequence of the second 0(1 /n) contri- 
bution is that while in the ideal case the probability 1/2 
is an upper bound, in the lossy case deviations from the 
leading term, 

p mal (i) = \ cxp(-2a; acot x) Y~2 ' ( 20 ) 

can be expected in cither direction. The functional form 
of Eq. ([20|) . plotted on Fig. [3l allows for a universal in- 
terpretation of the dependence of success probability on 
the transmission rate and the size of the search space 
through the combined variable x. For small losses we can 
use the approximation (fT7|) and conclude that the search 
efficiency depends only on the difference n/2 — e. The 
approximation is compared with the results of numerical 
calculations on Fig.[U We can observe the 0(l/n) accu- 
racy of the theoretical curves as expected, hence produc- 
ing poorer fits at smaller ranks. The positive deviations 
from the theoretical curves observable at low transmis- 
sion rates are due to the second 0(1/ n) term of Eq. (Q2 



V2 n ~ 1 acot(- I117/V2"- 1 ), 



(14) 



IV. DIRECTION DEPENDENT LOSS 



or the closest integer, as the time yielding the maximum 
probability of success. 



In the present section we no longer assume equal loss 
rates, and consider the schematically depicted loss model 



5 




FIG. 3: Probability of measuring the target state after the 
optimal number of iterations according to the approximation 
in Eq. I|20[) . The probability is plotted against the logarithm 
of x which is a combination of the rank of the hypercube n 
and the logarithmic transmission parameter e. 
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Rank of hypercube (n) 

FIG. 4: Maximum success probabilities for different ranks of 
hypercube (n) calculated using the theoretical approximation, 
and numerical simulations. The theoretical curves are drawn 
with continuous lines of different patterns, and the numeri- 
cal results are represented by points interconnected with the 
same line pattern and colour as the theoretical approximates 
corresponding to the same logarithmic transmission parame- 
ter e. 



on Fig. [5] with arbitrary rj d parameters. Because of the 
high symmetry of the hypercube graph, and the use of 
mainly identical multiports, we can neglect the position 
dependence of the transmission coefficients. The oper- 
ator T> describing the decoherence mechanism thus acts 
on a general term of the density operator as 

T>(\d,x) (d',x'\) = r\ d r\ d * \d,x) (d',x'\ + 8 xx >8 dd >r] d |0)(0| . 

(21) 

To describe the overall effect of this operator on a pure 
state, we re-introduce the linear decoherence operator in 



a more general form, 

D = Y j m\d)(d\®t i (22) 

d 

and use the notation {77} to denote the set of coefficients 
rj d . Due to the symmetry of the system, the sequential 
order of coefficients is irrelevant. With the re-defined 
operator the effect of decoherence reads 

V{q) = q' + (1 -Tr^) |0><0| , (23) 

where g = (ip\ is the initial state, and the non- 
vacuum part of the output state is g' — \ip') (ip'\, with 
= D\ip}- Therefore, we can again reduce our prob- 
lem to calculating the evolution of unnormalized pure 
states, just as in the uniform case, and use the non- 
unitary step operator U" = DU' with the more general 
noise operator. 

Telling how well the algorithm performs under these 
conditions is a complex task. First we give a lower bound 
on the probability of measuring the target node, based 
on generic assumptions. To begin, we separate the noise 
operator into two parts 

D = n + D', (24) 

where, for the moment, we leave < 77 < 1 undefined. 
As a consequence of Eq. (|2"2"|) the diagonal elements of D' 
arc [D'] dd = S d = rj d — rj, and the off-diagonal elements 
are zero. From Eq. ((23)) it follows that starting from a 
pure state \ipo), after t non-ideal steps the state of the 
system can be characterized by the unnormalized vector 
\ip'(t)), which is related to the state obtained from the 
same initial state by t ideal steps as 

W{t))=rfm)) + \r). (25) 
The expression of the residual vector |r) reads 
t 

|r) = W" 1 \m) • (26) 

fe=r 

To obtain the probability of measuring the target state 
I a; = 0) we have to evaluate the formula 

n-l 

Pn({v},t) = £ \rf(d, 0| m) +{d,0\ r)\ 2 . (27) 

Due to the symmetry of the graph and the coins, we use 
e.g. Eq. HU) and obtain (d, 0| ?/>(£)) w - sm(ui' t)/V2n. 
To obtain a lower bound on Pn({ij}^) we note that the 
sum is minimal if (d, 0| r) = const = K for every d (we 
consider a worst case scenario when all ( d, 1 r) are neg- 
ative) . Now we assume that the second term is a correc- 
tion with an absolute value smaller than that of the first 
term. For the upper bound on K, we use the inequality 

n-l 

Y,\(d,0\r}\ 2 <(r\r). (28) 

d=0 
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The norm of |r) can be bound using the eigenvalues of 
U, D, and D' . Let rj max = max {rj^d = 0, . . . ,n— 1} and 
<^max = max {\S e i\ \d ~ 0, . . . , n — 1}. Then we have 



k=l 



k-1 



V ??max ~ T] 



-fa 



t 

max 



»?*)• 

(29) 

Since U is unitary, its contribution to the above formula 
is trivial. Our upper bound on \K\ hence becomes |-fiT| < 

1 /v / "( 7 7max'5max/'7)('7max- 7 7*)/ (Vmax-v)- Combining the 
results, we obtain a lower bound on the probability for 
measuring the target node, 



p n m,t)>v 2t \^Ht) 



/max u max 



?? f/max - V 



n 



(30) 



where pn ^ (t) stands for the corresponding probability of 
the ideal (lossless) case. We maximize the lower bound 
with respect to the arbitrary parameter 77. The procedure 
can be carried out noting that 5 max = max{?j mra — 77, 77 — 
?7min}, thereby we find the maximum at 77 = fj = (rj max + 
??min)/2, yielding the formula 



Pn({v},t) > V 



21 



(imax/l) (?7maxA/) 



1 



(31) 

To interpret the formula (|3ip . we consider the two terms 
in the curly braces separately. The first term returns the 
success probability for uniform losses with transmission 
coefficient 77. The second term may be considered as a 
correction term that depends not only on some average 
value of the loss distribution, but also on its degree of 
non-uniformity in a way that is reminiscent of a mean 
square deviation. We observe that Eq. (|3"01) provides a 
useful lower bound only for {77} distributions violating 
uniformity to only a small degree. When the expression 
inside the curly braces becomes negative, the assumption 
made on the magnitude of the second term of Eq. (|2"T|) 
becomes invalid, and therefore the formula does not give 
a correct lower bound. 

The estimated lower bound (|3"0|) decreases with in- 
creasing degree of non-uniformity, in accordance with a 
naive expectation. However, as we shall show later, nu- 
merical simulations taking into account the full complex- 
ity of the problem provide evidence to the contrary: de- 
parture from uniformity can result in improved efficiency. 

Inspired by the appearance of the average loss rate in 
the lower bound (|3 1|) . we introduce the mean and the 
variance of the direction dependent losses, 



1 n-1 

(v) = -zZvd, and Q 



d=0 



-. n— 1 
d=0 



(32) 



By using the Taylor expansion of the success probability 
function p™ ax around the point r\d = (77) , the deviations 
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Average transmission ({»?)) 

FIG. 5: The coefficients for the (77) dependent second or- 
der term in the Taylor series expansion of jC ax ({?7}). The 
coefficients have been obtained by second order fitting to nu- 
merically generated values for a hypercube sized n = 8. The 
confidence of each fit is represented on the graph by an er- 
ror bar. It is clearly visible that the coefficients are always 
greater than 2~ n . Another feature, that is more suggestive 
on a linear scale, is that the points between two steps seem 
to align into straight lines, with their slopes increasing with 
(V)- 



from the uniform loss case can be well estimated at small 
degrees of non uniformity. Using the permutation sym- 
metry of p" lax we can express the Taylor series as 



max 



(M) 



max 



(( V )) + BQ 2 + CW 3 + 0(Sj), (33) 



where W 3 = ^/nj^k^l- We notice that Q may be re- 
garded as the mean deviation of {77} as a distribution, 
and hence it is a well-defined statistical property of the 
random noise. In other words, as long as a second order 
Taylor expansion gives an acceptable approximation, the 
probability of success depends only on the statistical av- 
erage and variance ((77), Q) of the noise and not on the 
specific values of {77}. Using numerical simulations, we 
have determined the values of B up to rank n = 10, and 
studied the impact of higher order terms. 

The second order Taylor coefficients were determined 
by fitting over the numerically obtained success proba- 
bilities at data points where the higher order moments of 
the loss distributions were small. An example plot of B 
is provided on Fig. [5J for a system n — 8. The higher or- 
der effects were suppressed by selecting the lowest values 
of W from several repeatedly generated random distribu- 
tions {77}. A general feature exhibited by all studied cases 
is that the second order coefficients satisfy the inequality 



B > 2" 



(34) 



It is remarkable that this tight lower bound depends only 
on the size of the system. The dependence of B on (77) is 
monotonous with discontinuities. We found the number 
of discontinuities to be proportional to the rank n. Our 
numerical studies have shown that the value of B before 
the first discontinuity is always a constant, and equal to 
the empirical lower bound ()34D . 
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FIG. 6: The relative improvement of maximum success prob- 
ability comparing direction dependent loss to uniform loss 
with identical average loss rates. The difference is measured 
as K ax (M) -Pn ax ({'?>)]Mr x ({?7}). All points were gener- 
ated at Q — 0.35 and the entire available domain for (rj). The 
values of higher moments of {77} are not restricted, therefore 
we see multiple points for certain ({rj) , Q) pairs. (Rank of the 
hypercube n = 9.) 

To plot the success probabilities corresponding to ar- 
bitrary random coefficients we used the pair of variables 
{i]) and Q. On these plots, the higher order terms cause 
a "spread" of the appearing curves. A sample plot is 
displayed on Fig. [5] where the relative improvement is 
compared to the uniform case, in percentages. We ob- 
serve a general increase of efficiency as compared to the 
uniform case with the same average loss rate. A general 
tendency is that for smaller values of (rj) the improve- 
ment is larger, interrupted, however, by discontinuities. 
These discontinuities closely follow those of the second 
order coefficient B. 

The numerical studies, involving the generation of 1000 
sets of uniformly randomly generated transmission coef- 
ficients for each of the systems of up to sizes n = 10, in- 
dicate that with the help of Eq. (|34|) the first two terms 
of the expansion Eq. can be used to obtain a general 
lower bound: 

pT x ({v}) > pT x ((v)) + 2-"Q 2 . (35) 

The inequality implies that the overall contribution from 
higher order terms is positive, or always balanced by the 
increase of B. The appeal of this lower bound is that it 
depends only on the size of the system N = 2™, and the 
elementary statistical properties of the noise ((t?), Q)- 
Therefore, together with the formula P0|) for uniform 
loss, a straight-forward estimation of success probability 
is possible before carrying out an experiment . 

Up to now, we concentrated on comparing the per- 
formance of the search algorithm suffering non-uniform 
losses with those suffering uniform loss with coefficient 
equal to the average of the non-uniform distribution. An- 
other physically interesting question is how attenuation 
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FIG. 7: The difference of the maximum success probabili- 
ties, in the presence of direction dependent loss with coeffi- 
cients {r/}, and in the presence of uniform loss with coeffi- 
cient 77 max = max{?7}. Physically, the non-uniform case can 
be obtained from the uniform case by introducing proper at- 
tenuation. We use the second moment Q as a measure of 
deviation from the non-attenuated case. The vertical axis 
shows [Pn ax ({rj}) -p" lax (r?max)]/p™ ax (?)max) as a percentage. 
We observe a systematic improvement for higher values of Q. 
The plot corresponds to rank n = 7 and rymax = 0.996±0.001. 



alone affects search efficiency. We can formulate this 
question using the notations above as follows. Consider 
a randomly generated distribution {77} and compare the 
corresponding success probability with the one generated 
by a uniform distribution with transmission coefficient 
Vmax = max{?7}. We chose Q as a measure of how much 
an ?7max uniform distribution needs to be altered to ob- 
tain {77}, and made the comparisons using the same set 
of samples. A typical plot is presented on Fig. [7J It 
appears that as we start deviating from the original uni- 
form distribution, an initial drop of efficiency is followed 
by a region where improvement shows some systematic 
increase. However, it is still an open question, whether it 
is really a general feature that for some values of Q the 
efficiency is always increased. On the other hand these 
plots provide clear evidence that for a significant number 
of cases the difference p™ ax ({^}) — Kf ax (?7max) is posi- 
tive. In other words, rather counter-intuitively, we can 
observe examples where increased losses result in the im- 
provement of search efficiency. Since the time evolution 
with losses is non-unitary, the improvement cannot be 
trivially attributed to the fact that the Grover operator 
is not the optimal choice for the marked coin. 



V. PHASE ERRORS 

In the present section we discuss another type of errors 
typically arising in optical multiport networks. These er- 
rors are due to stochastic changes of the optical path 
lengths relative to what is designated, and manifest 
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as undesired random phase shifts. Depending on how 
rapidly the phases change, we may work in two comple- 
mentary regimes. In the "phase fluctuation" regime the 
phases at each iteration are different. These errors can 
typically be caused by thermal noise. In the "static phase 
errors" regime, the undesired phases have slow drift such 
that on the time scale of an entire run of the quantum 
algorithm their change is insignificant. The origin of such 
errors can be optical clement imperfections, optical mis- 
alignments, or a slow stochastic drift in one of the ex- 
perimental parameters. Phase errors in the fluctuation 
regime have been studied in Ref. [22j for walks on N 
dimensional lattices employing the generalized Grover or 
Fourier coin. The impact of a different type of static error 
on the SKW algorithm has been analyzed in Ref. [25| . 

To begin the formal treatment, let F denote the oper- 
ator introducing the phase shifts, and write it as 



</../• 



This operator is unitary, hence the step operator 
U(M) = SF({<p})C, 



(36) 



(37) 



that depends on the phases {<pdx\d = 0..n— I, x = 0..2" — 
1} is unitary as well. In case of phase fluctuations, at 

each iteration t we have the parameters ip^ such that all 

rdi are independent random variables for every d, x and 
t, according to some probability distribution. In case of 
static phase errors, ip^. and tp^J are considered to be 
the same random variables for every pair of d and x. 

The formalism of Ref. [22j can be applied to the walk 
on the hypercube, and extended to the case of non- 
uniform coins and position dependent phases. Namely, 
using the shorthand notations D = {0, 1, 2, . . . , n — 1} 
and 



E(fc,l)=0e.„ 

3=1 

the state after t iterations can be expressed as 



(38) 



WM,*)> 



f_i\S*t*o e itp(a,x ) 

xZ Xt (a,x )\ai,x ®E{t,l)), (39) 



where 



S Xt (a, x Q ) 
t-i 



a j + 1 S Xt ®x ,E(t,j+l) j , 

3 = 1 

(40) 



and <p(a,x ) = Yfj=i ( Pa + ,x iE(t,j-i)- For thc standard 
SKW algorithm, the coin matrices are C^J, = 2/n — 8 aa i 



and C^J, = —6 aa ', however, the SKW algorithm is re- 
ported to work with more general choices of operators 
Co/i 1- 

For the following study, we express the probability 
of finding the walker at position x after t iterations as 
the sump n (x, {</?},*)= pi(x, {cp}, t) +p„(x, {ip}, t), such 
that the incoherent and coherent contributions are 



aeD* 



(41) 



^ E ^a^xMr^M^m 



where x t = Xt ffi x. The appearing phase factors are 
and S St (a) = E Xt (a, E(t,l)), i.e. 

3 = 1 

(44) 

Note, that when the probability of finding the walker at 
the target node xt is to be calculated we must set x = Xt , 
therefore, we have x t — 0. 

In the following we shall show that the incoherent con- 
tribution is constant. 



(45) 



for any two unitary coins Co/i- Consequently, p 1 ^ is con- 
stant also for balanced coins such as those in Eq. ([5]). The 
summations in Eq. (|4ip can be rearranged in increasing 
order of indices of cij , yielding 

p I n (x,{<p},t) = 

71—1 

— V lc (0) + f^ (1) - ^ (0) l 6- 



n2 



n-1 



a t _i=0 
n-1 



E ^ iat + [C( 1 )-C(°)] at _ iat <5 £t , i;(t „ 1 , 1 



<U=0 



(46) 



Since E(t — 1,1) depends on a,j only when j < t — 1, 
and due to the unitarity of the coins (a t _i |CgCo|at_i) = 
(at_i|C|Ci|at_i) = 1, the summation over at can 
be evaluated and we obtain 1. Hence, we see that 
p n {x, {(p}, t) = p n (x, {(p},t— 1), and this implies Eq. (|4"5"|) 
by induction. 

The average probability p n (x,t) of finding the walker 
at node x is obtained by averaging the random phases 
according to their appropriate probability distribution. 
Using Eq. (|4"5|) this probability can be expressed as 

Pn(x, t) = ( Pn (x, M, t)) = — + M, i)) , (47) 
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Number of steps (/) 

FIG. 8: The averaged (1000 samples) probability of measuring 
the target node, when n = 6 and Aip = 3°, 6°, 9°, 12°. The 
tendency of the success probability to a constant, non-zero 
value can be observed on this numerically obtained plot. It 
is also observable that a larger variance results in a smaller 
asymptotic value. 



0.4 




Number of steps (t) 

FIG. 9: Time dependence of success probabilities for two dif- 
ferent phase configurations, numerically calculated for a sys- 
tem of rank n — 6. The difference in frequencies of the major 
oscillations is clearly observable for larger times. 



where (...) denotes taking the average for each random 

variable ip^l m case 01 phase fluctuations, and for each 
ifdx in case of static phase errors. It is reasonable to as- 
sume that each random variable has the same probability 
distribution. To analyze the impact of phase errors on the 
search efficiency, we study the behaviour of the coherent 
term (p%(x, {<p}, t)) for different random distributions. 

In case of phase fluctuations characterized by a uniform 
distribution, the coherent term immediately vanishes and 
we obtain p n (x,t) = 1/2™. This case can be considered 
as the classical limit of the quantum walk. Therefore, we 
conclude that the classical limit of the SKW algorithm is 
not a search algorithm, independently of the two unitary 
coins used. 

Assuming a Gaussian distribution of random phases 
is motivated by the relation of each phase variable p> to 
the optical path length. The changes in the optical path 
lengths which introduce phase shifts are not restricted to 
a 2tt interval. In what follows, we assume that the ran- 
dom phases have a zero centered Gaussian distribution 
with a variance Aip. 

We arrive at the classical limit even when the phase 
fluctuations have a finite width Gaussian distribution, 
simply by repeatedly applying the time evolution oper- 
ator U{ip}. For such Gaussian distribution, the coher- 
ent term exhibits exponential decrease with time, a be- 
haviour also confirmed by our numerical calculations. 

In the static phase error regime the mechanism of can- 
cellation of phases is different than in the fluctuation 
regime, and more difficult to study analytically. For uni- 
form random distribution we expect a sub-exponential 
decay of the coherent term to zero. For a zero centered 
Gaussian distribution with variance Aip we performed 
numerical simulations using the standard two coins of 
Eq. ©. 



The numerical results for the success probability 
Pn{xt,t) for several values of A<p are plotted on Fig. [5] 
The data points were obtained by calculating success 
probabilities for 1000 randomly generated phase config- 
urations and taking their averages at each time step t. 

By studying the repetition of the random phase config- 
uration we come to several remarkable conclusions. First, 
the time evolution of the success probability tends (on 
a long time scale, t > tj) to a finite, non-zero constant 
value. Consequently, being subject to static phase errors, 
the SKW algorithm retains its utility as search algorithm. 
Second, the early steps of the time evolution are charac- 
terized by damped oscillations reminding of a collapse. 
Third, the smaller the phase noise the larger is the long 
time stationary value to which the system evolves. We 
have plotted the stationary values obtained by numer- 
ical calculations, against the rank of the hypercube on 

Fig.nni 

Better insight into the above features can be gained by 
examining the shape of the individual runs of the algo- 
rithm with the given random phase configurations. As it 
can be seen on Fig. [51 the success probabilities for dif- 
ferent runs display the typical oscillations around a non- 
zero value. They differ slightly in their frequencies de- 
pending on the random phases chosen, hence when these 
oscillations are summed up we get the typical collapse be- 
haviour. Also, since these frequencies continuously fill up 
a band specified by the width of the Gaussian, we expect 
no revivals to happen later. For higher order hypercubes 
the success probability drops almost to zero for already 
very moderate phase errors, resembling a behaviour seen 
on Fig. H 
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FIG. 10: The long time stationary values of the success prob- 
ability (obtained by averaging over 1000 samples) against the 
size of the search space, for A<p = 0°, 3°, 6°, 15°. 



VI. CONCLUSIONS 

We studied the SQRW implementation of the SKW 
search algorithm and analyzed the influence on its per- 
formance the two most common type of disturbances, 
namely photon losses and phase errors. Our main result 
for the photon loss affected SQRW search algorithm is 
that the introduction of non-uniform distribution of the 



loss can significantly improve the search efficiency com- 
pared to uniform loss with the same average. In many 
cases, even the sole increase of losses in certain directions 
may improve the search efficiency. Mostly based on nu- 
merical evidence we have set a lower bound for the search 
probability as a function of the average and variance of 
the randomly distributed direction dependent loss. 

We concentrated our analysis on two complementary 
regimes of phase errors. When the system is subject to 
rapid phase fluctuations, the classical limit of the quan- 
tum walk is approached. We have shown that in this 
limit the SKW algorithm loses its applicability to the 
search problem for any pair of unitary coins. On the 
other hand, we showed that when the phases are kept 
constant during each run of the search, the success rate 
does not drop to zero, but approaches a finite value. The 
effect in its mechanism is reminiscent to exponential lo- 
calization found in optical networks [30]. Therefore, in 
the long-time limit, static phase errors are less destruc- 
tive than rapidly fluctuating phase errors. 
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